#include "solver_stable.h"

#include <math.h>
#include <memory.h>

#include "FuncSolver.h"

double* solve_stable(double* u0, Area area, double epsilon, TwoDimFunc f)
{
	task_scheduler_init init;

    int n = area.n;
    double a = area.a, b = area.b;
    double* u_cur = new double[n * n];
    double* u_prev = new double[n * n];
    double h = (b - a) / n;
    double error;
    double start, finish;

    double* f_cur = new double[n * n];

    start = PortableGetTime();
	/*
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            f_cur[i * n + j] = f(a + i * h, a + j * h);
        }
    }*/

	parallel_for( blocked_range<int>(0,n,1), FuncSolver( f_cur, f, n, a, b, h ) );

    finish = PortableGetTime();
    printf("Compute function time: %.4f\n", finish - start);

    memcpy(u_prev, u0, n * n * sizeof(double));
    memcpy(u_cur, u0, n * n * sizeof(double));

    start = PortableGetTime();

	tick_count st,fh;

	double f1 = 0,
		f2 = 0;
    do
    { 
			st = tick_count::now();
		parallel_for( blocked_range2d<int>(1,n-1,1,n-1), EquSolver( u_cur, u_prev, f_cur, n, h ) );
			fh = tick_count::now();
			f1 += (fh-st).seconds();

			st = tick_count::now();
		error = abs(u_cur[0] - u_prev[0]);
		ReduceErr r( error, u_cur, u_prev );
		parallel_reduce( blocked_range<int>(1, n*n - 1), r );
		error = r.Result();
			fh = tick_count::now();
			f2 += (fh-st).seconds();

        double* tmp;
        tmp = u_cur; u_cur = u_prev; u_prev = tmp;

    } while (error > epsilon);

    finish = PortableGetTime();

	printf("----------Time of solver: %.4f\n", f1);
	printf("----------Time of checker: %.4f\n", f2);
    printf("Solve time: %.4f\n", finish - start);

    delete[] u_cur;
    return u_prev;
}